Heat transport and thermal rectification in quasi-one-dimensional systems. 
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In this work we investigate heat conduction along a ladder-model conformed by two coupled one 
dimensional lattices with different anharmonicity. We study how the interchain coupling modifies 
the thermal properties of the isolated systems. For a large enough coupling strength, we demon- 
strate that a harmonic lattice interacting with an anharmonic one is able to support a linear thermal 
gradient when it is connected to two heat reservoirs at different temperatures. We estimate this 
critical coupling by applying the self-consistent phonon theory (SCPT) to the anharmonic counter- 
part. By exchanging the heat baths connections between the harmonic and the anharmonic chains, 
our results show that the coupled system reveals as a thermal rectifier. 
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I. INTRODUCTION 

The current development of nanomaterials for molecu- 
lar electronics has triggerd the study of heat conduction 
in lovi'-dimensional systems from new perspectives. It is 
well known that the derivation of the thermal proper- 
ties in one and two dimensional systems is a controver- 
sial topic which has been deeply studied in the last two 
decades >ii^ In the case of perfect harmonic lattices, the 
lack of phonon-phonon interactions leads to no thermal 
resistivity, which gives rise to an infinite thermal conduc- 
tivity in the thermodynamic limit. This means, Fourier's 
law is not valid in these systems and heat conduction ex- 
hibits anomalous behavior. Nonintegrability and an ex- 
ternal substrate potential constitute the sufficient condi- 
tions to completely change this scenario.'^ In such a case, 
the conservation law of total momentum is violated and 
Fourier's law reveals as valid. Recently, the most studied 
proposals are the Frenkel-Kontorova (FK))^"— and the cj)'^ 
potentialsj^"— Several authors have revisited the unusual 
thermal properties of one dimensional systems in order 
to propose different applications such as thermal recti- 
fiers or thermal transistors J^""— Some of these studies 
have been extended to two and three dimensional struc- 
tureai^iii and complex networks^®' where it has been 
addressed the relevant influence of the coupling among 
several one dimensional lattices. 

More interestingly, the analysis of the mechanisms me- 
diating energy flow in low dimensional biomolecules is a 
fundamental issue for the understanding of many bio- 
logically relevant functions. The electronic and vibra- 
tional degrees of freedom of biomolecular systems, spe- 
cially those containing helix structures , i.e. double-helix 
DNA or a-helices in proteins, can be modeled as lad- 
der structures of coupled one-dimensional lattices<^"— 
Recently, the thermal conductivity of double-stranded 
molecules has been studied in Ref. [13 where the inter- 
chain interaction between two identical chains has re- 
vealed as a positive or negative effect on heat conduction, 
depending on the strength of the nonlinearity present in 
the system. 



In this paper we address some interesting issues about 
heat transport along double-stranded molecular systems 
which remain open. In particular, we consider a ladder- 
model where two different lattices, a harmonic and an 
anharmonic one, are coupled by harmonic forces. In view 
of the clearly different thermal properties of these two 
subsystems when they are isolated, our work analyzes 
how the coupled system behavior is affected depending on 
the coupling strength. In addition, due to the asymmetry 
present in the system, we address the possibility of heat 
rectification. 



II. MODEL 

In this section we will present the theoretical formal- 
ism we will be dealing with in our study on heat trans- 
port along double-stranded molecules. In particular, we 
consider a system conformed by a harmonic one dimen- 
sional lattice which is coupled by harmonic forces to an 
anharmonic one. The dimensionless Hamiltonian of such 
a system is written as follows: 
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T-La and Hh describe the dynamics of the subchains, 
the harmonic (H) and the anharmonic (A) one respec- 
tively, each one consisting of n ~ 1 . . . N sites. Both 
are affected by a harmonic potential W{xn,Xn-i) — 
^{xn — Xn-i)"^ and chain (A) is considered within the 
Frenkel-Kontorova model. Thus, chain (A) is affected 
by an on-site potential V{xn) = ^^cos27rx„, whose 



2 



strength is defined by Vq . For the sake of simphcity, we 
will consider the same elastic constant for both chains, k, 
in terms of which the harmonic interchain coupling kint 
is expressed. The latter interaction is considered within 
the Hamiltonian T-Lmt , whose characteristic oscillator fre- 
quency is ujQ = k/m . Notice that hereafter the period 
of the external potential, b, is taken as the length scale 
of the system and the magnitude kb^ will be the energy 
unit 4 

Fig. [1] presents a schematic view of the system under 
study. 
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FIG. 1: Schematic view of a fragment of the ladder-model 
under consideration. The harmonic lattice is connected at 
left and right edges to two independent heat reservoirs at 
temperatures Tc and Th respectively. 

Our aim is to study the formation of a thermal gra- 
dient along this double-stranded structure by connecting 
two independent heat reservoirs at both ends of the har- 
monic lattice. Heat baths effects will be included in the 
dynamics of the edge sites of the system, yi and un, hy 
using stochastic Langevin equations.— This will modify 
the Newton's equations of motion of the affected sites, 
n = 1 and n = A^, as follows; 

^ = -W'{yn,yn-i) ~W'{yn,yn+i) - 1^ + fn{t) , 

(2) 

where 7 is the coupling between the system and the 
stochastic bath and the prime indicates derivative with 
respect to Random forces, which include temperature 
effects, are defined as: 



(3) 



The relationship between this dimensionless tempera- 
ture T and the real one is given by T = KbT^ /kb"^, 
where Kb refers to the Boltzmann constanti^ The tem- 
perature of the baths T will be different for the left cold 
{Tc) and right hot {Th) sides of the chain. 

Exploiting the equations of motion, the local heat fiux 
is defined by the continuity equation as: 
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Molecular dynamics simulations are performed by the 
numerical method proposed by Greenside and Helfand 
as a correction to the Runge-Kutta method for stochas- 
tic equations (3o4s2g),^^"^^ by considering a long enough 
integration time such that the stationary state is estab- 
lished. The time step is 5t = 10^^ and the friction con- 
stant of the baths is set to 7 = 0.5 in all simulations. 
In the final state the time averaged heat flux reaches a 
constant value along the system such that J = (Ji) = 
... = (Jat). Thus, the thermal conductivity for a finite 
system can be calculated as k = JN/{Tc — Th). Notice 
that K will be size-independent in the case of normal heat 
transport. Similarly, in the steady state the time aver- 
aged temperature calculated as T„ = (i^-l-y^) will reach 
the stationary thermal profile. 



III. SINGLE ANHARMONIC CHAIN 

To better understand the infiuence of the coupling to 
nonlinear degrees of freedom, first we would like to per- 
form a preliminary study of a single Frenkel-Kontorova 
chain connected to two independent heat reservoirs at 
temperatures Tc — Tm — ST and right Th = Tm + ST at 
the left and right edges of the chain. Thus, we integrate 
the Newton's equations derived from "Ha and heat baths 
effects are included in the dynamics of sites xi and xn 
according to Eq. [21 

We have performed molecular dynamics simulations for 
different values of Tm, ST = 0.05 and Vq = 8. Once the 
steady state is reached, the stationary local temperature 
profile is calculated as T„ = (i^) which is shown in left 
panels of Fig. [2l 

As it is already known, our results in Fig. [21(b) and (c) 
show that if the FK-potential strength is large enough re- 
spect to the temperature Tm, a well-defined thermal gra- 
dient is formed. Furthermore, in such a case, the Fourier 
heat law is valid and the total heat flux along the system 
does not diverge in the thermodynamic limit contrary to 
the situation occurring in harmonic systems. Fig. [H^a) 
demonstrates on the other hand that Tm has to be high 
enough to be able to activate the vibrational states of 
the anharmonic band. This can be understood by invok- 
ing a Landauer-like equation to calculate heat flux in the 
ballistic regimej^ 



J 



huj[nc{i^) — nH{td)]T{uj)duj. 
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Here t{ll)) refers to the transmission coefficient, being 
uj the frequency of every vibrational state, and the 
mode distribution of the heat baths are considered clas- 
sical within the Maxwell-Boltzmann statistics nii^L{uj) = 
exTp{~huj / K bT^ j^) . In the limit of linear response, 
ST << Tm, the thermal conductivity can be written asi^ 



K = 2^ ^ J exp{-x)T{x)dx. 



(6) 
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FIG. 2: Left panels show the stationary temperature profile 
calculated by numerical simulations in a FK lattice of A'^ = 50 
sites connected to two independent heat reservoirs at different 
temperatures, Th,c ~ Tm ± 0.05. Right panels present the 
comparison between the numerical DOS calculated in equi- 
librium with the theoretical frequency band obtained accord- 
ing to the SCPT (shadowed region). The solid line accounts 
for the weight function F{Tm)- The temperatures considered 
in the simulations are a) Tm = 0.10, b) Tm = 0.15 and c) 
Tm = 0.20. The inset shows the frequency band of states 
calculated by SCPT as a function of temperature and Vb = 8 
in the shadowed region. 

The integral in Eq. [HI being x = huj / K bT^j , is given 
by the product between the transmission coefficient and 
a weight function F(Tm) = a;^exp(— a;), which accounts 
for bath thermal effects. This means that F{Tm) de- 
fines which vibrational system modes can be excited by 
the heat reservoirs and contribute to the heat fiux. Right 
panels of Fig. [5] show the numerical density of vibrational 
states (DOS), calculated by Fourier transforming the ve- 
locity autocorrelation function of equilibrium molecular 
dynamics simulations, namely, Tc = Th — T^/. We also 
plot the weight function F{Tm) for the Tm considered in 
the simulations. It is clear that for low Tm there is very 
few states under the curve F{Tm), and therefore, only 
those states are involved in heat transfer along the sys- 
tem, see left panel of Fig.l^Ja). The lack of heat carrying 
modes gives rise to a stationary thermal profile deviated 
from the linear thermal gradient and thus, Fourier's law 
is not expected to be valid. 

To get further insight into the behavior of the nonlinear 
system the self-consistent phonon theory (SCPT) can be 
considered. This approach consists in replacing the an- 
harmonic potential for a harmonic approximation such 
that the new effective harmonic strength is temperature- 
dependent as follows >^i2& 

- ^ cos(2^x„) ^ , (7) 

By performing a variational study,^^ the effective in- 



teraction U{T) for the Hamiltonian Ha can be obtained 
by solving the following transcendental equation: 

U{T) = Vo exp 8 

Once U (T) is known, the dispersion relation of the sys- 
tem, w\{p,T) = U{T) + 2(1 - cos{p)), being p e [0,7r], 
and the DOS for the effective harmonic system can be 
defined analytically. Hereafter u refers to frequencies ex- 
pressed in units of the fundamental frequency ujq = k/m 
of the system. The frequency band of the SCPT for 
Vb = 8 as a function of the temperature is presented 
in the shadowed region of the inset of Fig. [2] The com- 
parison between the frequency band defined by {p, T) 
and the numerical DOS calculated in equilibrium, is also 
shown in the right panel of Fig. [3] by considering the so- 
lution of Eq. [8] at T = Tm. 

Figure [5] shows from top to bottom that the anhar- 
monicity shifts the DOS to lower frequencies when Tm is 
increased, as predicted by the SCPT approach. However, 
it is well known that the SCPT approximation is not able 
to predict the broadening of the frequency band when 
thermal effects become more relevant at higher tempera- 
tures, see Fig. [SJc). In view of these arguments, we want 
to stress two conclusions which will be taken into account 
hereafter. On one hand it is clear that T^v/ in the system 
should be high enough to activate the anharmonic modes 
in order to create a well formed thermal gradient along 
the lattice in presence of the thermal baths. However, 
in order to have an accurate description of the thermal 
properties based on the SCPT, Tm should be low enough 
so that thermal effects are kept within this approach va- 
lidity. For the sake of clarity we will focus on temperature 
regimes which fulfilled these two requirements hereafter. 

IV. LADDER-MODEL 

Now that we have reviewed some of the characteristics 
of an anharmonic one dimensional lattice, we focus on the 
main system of interest in this work, the ladder-model 
presented in Fig. [U whose Hamiltonian is described by 
Eq. [1] Notice that the heat reservoirs at temperatures 
Tc and Th are connected to both edges of the harmonic 
chain. Therefore, the dynamics of sites yi and j/jv will be 
described by Eqs. [H accordingly. 

As previously mentioned, a well formed thermal gra- 
dient cannot arise along an isolated harmonic chain con- 
nected to two heat baths at different temperatures. How- 
ever, its thermal properties are expected to change by 
switching on its interaction with the anharmonic chain 
in our system. Thus, by increasing kint the effect of 
the anharmonicity will be intensified in the whole ladder- 
system. Our results demonstrate that for a large enough 
interchain interaction, the vibrational spectra of both 
chains become mixed, and a linear thermal profile will 
arise in the system. This transition is shown in left pan- 
els of Fig. [3] for a system of iV = 100 sites, where the 
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anharmonicity is Vq = 5 and two interchain interactions, 
kint = 0-05 and hint = 1-00 are considered. Figure [3fc) 
also shows that the thermal conductivity along the lattice 
decays as a function of kint- This is a result of the fact 
that the contribution from the anharmonic chain to the 
double-lattice dynamics is stronger for a larger coupling 
and therefore, the thermal resistivity is higher. More 



bands defined by fci„t) and uj\{p,Tm) will have 

spectral overlap. According to several simulations, we 
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FIG. 3: Left panels show the stationary temperature profile 
calculated by numerical simulations in a ladder- model (Fig.H]) 
of N — 100 sites connected to two independent heat reservoirs 
at different temperatures, Tc = 0.1 and Th = 0.2. The tem- 
perature profile is plotted for the harmonic T„{H) and the 
anharmonic Tn{A) chain, and the coupled system Tn. The 
considered parameters are Vo = 5 and a) ki„t ~ 0.05 and b) 
kint = 1.00. Right panels show the finite-size thermal conduc- 
tivity as a function of c) the coupling ki„t and d) the system 
size A'^ for kint = 0.05 and kint = 1.00. 

interestingly, the presence of the thermal gradient for a 
large enough coupling is associated to a finite-size ther- 
mal conductivity which does not depend on the system 
size, and therefore, Fourier's law is expected to be valid, 
see Fig. M^d). 

In order to give an estimation of the threshold inter- 
action necessary to create a linear thermal gradient 
in the ladder-model, we will consider a fully harmonic 
ladder- system, whose vibrational spectrum consists of 
two bands, an acoustic and an optical one. In partic- 
ular the optical band will be at fci„t) — 2kint + 
2(1 — cos(p)), namely, it will be shifted to higher frequen- 
cies by increasing the kint- Note that in the real system 
under consideration (Fig. [T]) for small couplings the an- 
harmonic band will play the role of the optical band in 
the harmonic ladder-model. Once the interaction is large 
enough so that the anharmonic lattice can be affected by 
the heat reservoirs and thus, it can contribute to establish 
the thermal gradient it will present a band of anharmonic 
states at frequencies u}\{p,Tm) — U{Tm) + 2{1~cos{p)). 
This means that, the higher the temperature affecting 
the anharmonic lattice is, the lower the frequencies of 
the activated anharmonic modes are, see Sec. lIIIl There- 
fore, we will define the threshold coupling such that the 
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FIG. 4: Left panels show the stationary temperature profile 
for the harmonic Tn(H) and the anharmonic chain Tn{A) cal- 
culated by numerical simulations in a ladder-model (Fig. [T|) 
of = 100 sites connected to two heat reservoirs at Tc =0.1 
and Th = 0.2. Right panels present the theoretical DOS of a 
single FK chain affected by a mean temperature Ta accord- 
ing to the SOFT, as well as the DOS of the optical band of a 
fully harmonic ladder-model for several kint- The considered 
parameters are Vo = 8 and a) kint ~ 0.1, b) kint = 0.25 and 
c) kint = 1.75. 

have concluded that the overlap condition between the 
bands 0J%(p,kint) and u}\{p,Tm) depends on the tem- 
perature present in the system. This is related to the 
thermal broadening effects mentioned in Sec. IIIII When 
these effects, not included in SOFT, are negligible, we 
have numerically tested that the thermal gradient arises 
if approximately half of the states are common to the two 
considered bands. Thus, the interaction threshold can be 
estimated by: 



k* — 



U{Tm) - 2 



(9) 



Figure |4] shows the formation of the thermal gradient 
along the ladder-model by increasing kint for a system 
of iV = 100 sites and Vq = 8. The right panels present 
the theoretical DOS corresponding to the band oj\{p, Ta) 
of a single anharmonic lattice affected by an averaged 
temperature Ta = obtained in the simu- 

lations for each coupling. Additionally, the DOS of the 
theoretical optical band for the fully harmonic ladder- 
model uijjip, kint) for the considered couplings is plotted. 
In Fig. [SJc) , it is shown a well- formed thermal gradient 
along the system for the coupling k*^^ = 1.75, for which 
the DOS of ijjjip, fc*„J and u}\{p, Ta) overlap in 50% of 
their area. In Fig. [5] we show a good agreement between 
our prediction of fc*„( according to Eq. IH] and the thresh- 
old interaction obtained by numerical simulations for a 
system size of = 100 sites and different values of Vq. 
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FIG. 5: Critical coupling necessary to create a thermal 
gradient in the ladder-system (Fig. \Q of N = 100 sites as a 
function of the anharmonicity Vq. Solid line represents the 
theoretical prediction Eq.[5]in comparison with the values for 
k*„t numerically obtained, solid dots. In all cases Tc = 0.1 
and Th ~ 0.2 have been considered 
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V. HEAT RECTIFICATION 

In previous works it has been well established that 
symmetry breaking and nonlinearity are the necessary 
ingredients in a system to support heat rectification J^"— 
In this regard, we propose a new thermal rectifying mech- 
anism present in the ladder-model under consideration. 
On the one hand, it is clear that nonlinearity comes from 
the anharmonicity of one of the chains. On the other 
hand, the symmetry breaking occurs by considering dif- 
ferent connections for the heat reservoirs, in a similar way 
which has been considered in the case of charge transport 
in DNA.^^ Figure ini shows the two considered baths con- 
figuration where Tjj and Tc refer to the temperature of 
the hot and the cold heat reservoirs as previously. Notice 
that in Fig. |6|a) the hot reservoir is connected to the 
left edge of the anharmonic chain while in Fig. ^h) the 
position of the baths connection is reversed. This means 
that the heat current flows in opposite directions for both 
cases but more interestingly, our results show that these 
currents reach different values in the stationary state, and 
therefore the system behaves as a thermal rectifier. 

We perform molecular dynamics simulations by con- 
necting the ladder-model described by the Hamiltonian 
Eq. [1] to two heat reservoirs. In Fig. [6ja) the dynam- 
ics of sites xi and ijn will be described by Eqs.H] under 
Langevin baths at temperatures Th = 0.2 and Tc = 0.1 
respectively. In Fig. ^h) the connections are reversed 
and sites j/i and xjv are the one affected by Eqs. [2] 
accordingly. In order to analysis the main features of 
the heat rectification, we will numerically calculate the 
thermal conductivity, a) = J+N/{Th — Tc) and b) 
K_ = J^N/{Tc — Th), for several couplings kinf No- 
tice that the interchain interaction has revealed as the 
key parameter to change the thermal properties of the 
ladder-model in Sec. [31 These results are presented in 



FIG. 6: Schematic view of a fragment of the ladder-model 
under consideration connected to two independent heat reser- 
voirs at Th and Tc in two different configurations. Notice 
that in a) the hot bath is connected to the anharmonic chain 
and in b) this connection is made to the harmonic one. Ac- 
cordingly heat current flows in opposite directions, J+ and 
J_, respectively. 



Figure EJa) for Vq = 5 and two different system sizes 
N = 100 and N = 200. 

Our results establish the existence of heat rectification, 
> K_, for a certain range of the interchain inter- 
action kint- The difference between both conductivities 
A — Kjf. — K- reaches a maximum A^ax for a particular 
coupling fc°jj 0.25, where we find a strong rectification 
effect of 35%. Qualitatively we can understand this recti- 
fying behavior by the following handwaving argument. It 
is clear that if there is no coupling between both chains, 
the system cannot reach any stationary state, since the 
heat current is not able to flow from one bath to the 
other. However, when the coupling is switched on, new 
channels connecting both chains allow for the heat flow 
from the hot reservoir to the cold one and the thermal 
conductivity increases. Due to the coupling in a) the 
heat current flows from a bad conductor, the anharmonic 
chain at Th, to a good one, the harmonic chain at Tc, 
while in b) the situation is the opposite one, see Fig. |6l 
Therefore, in both cases k+ and k_ increases with k^t 
but since the heat current is deviated to a more efficient 
heat conductor in a), is larger than k_. In order 
to support this argument Fig. EJb) shows that A^ax in- 
creases with the system size TV, this means with the num- 
ber of coupling channels, as expected. 

Once a good connection between both baths is estab- 
lished, if kint > k^nt according to Sec. [3l the effect of the 
anharmonicity will be intensifled in the ladder-system. 
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FIG. 7: a) Thermal conductivity for the two cases presented 
in Fig. [6l k+ and function of the interchain couphng 

with Vo ~ 5 and for two different system sizes A'^ = 100 and 
A'' = 200. b) System size dependence of the thermal con- 
ductivity difference Amax for Vo = 5. c) Normalized thermal 
conductivity difference A/ Amax as a function of the coupling 
kint for TV = 100 and two different anharmonicities Vo ~ 5 
and Vo = 8. 



and therefore the conductivity is expected to decrease, 
see Fig. ISjc). When the coupling is so large, that the 
vibrational spectra of both chains are completely mixed, 
the rectifying effect is expected to disappear and thus 
~ In Fig.[6l^a) we shovif that this situation occurs 
for a fc-, ^ 1.0. 

Last, we would like to address, how the anharmonicity 
Vq affects the main features of heat rectification in the 
system. For the sake of comparison Fig. IHl^c) shows the 
magnitude A/Aj^ax as a function of fcj„( for two differ- 
ent anharmonicities, Vq = 5 and Vq — 8. Our results 
show that the qualitative rectification behavior does not 
depend on Vq but its main features are shifted to larger 
couplings kint- In a similar way, in Sec. |3] it was al- 
ready demonstrated that the threshold interaction 
for which the vibrational spectra of the subchains are 
mixed in the system, increases with the anharmonicity 
Vq, see Fig. [5] It is to be mentioned that it is not trivial 



to establish a direct relationship between k!^^ and fc*„j 
since the baths connections and thus, their thermal ef- 
fects are different in Fig. [T] and Fig. [51 However, accord- 
ing to our simulations it turns out that k"^^ > fc*„t, and 



therefore, we can predict that the main heat rectification 
effects will appear for k^t < k*^^. 



VI. CONCLUSIONS 

In this paper, we have studied heat conduction along 
a ladder-model consisting of a harmonic and an anhar- 
monic one dimensional chain coupled by harmonic inter- 
actions. We have analyzed how the thermal properties 
of the coupled ladder-system depends on the strength of 
the coupling. In particular, by connecting two indepen- 
dent heat reservoirs to the harmonic counterpart of the 
system, we have demonstrated that for a large enough 
interchain interaction kint, a thermal gradient can be 
formed along the system and the thermal conductivity 
will remain constant for increasing size systems. This 
means that the ladder-system will present normal heat 
transport contrary to the case of an isolated harmonic 
system. We have estimated the threshold interaction of 
this transition by considering the self-consistent phonon 
theory for a single anharmonic chain in comparison to a 
fully harmonic ladder-model. Our estimations has been 
proven to be in good agreement with numerical results 
based on molecular dynamics simulations. 

Taking advantage of the nonlinear effects present in our 
system, it was shown that heat current along the ladder- 
model behaves differently if the heat baths connections 
are interchanged between the harmonic and the anhar- 
monic chain as in Fig. [S] This mechanism reveals strong 
heat rectification effects of more than 30%. The main 
qualitative features of the rectifying device do not depend 
on the system size nor on the anharmonicity present in 
the system. However, the maximum rectification rates 
increases with the system size and shifts to larger cou- 
plings in the case of stronger anharmonic effects. 
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